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Abstract. The focus of the present study is the modified Buckley-Leverett 
(MBL) equation describing two-phase flow in porous media. The MBL equa- 
tion differs from the classical Buckley-Leverett (BL) equation by including a 
balanced diffusive-dispersive combination. The dispersive term is a third order 
mixed derivatives term, which models the dynamic effects in the pressure dif- 
ference between the two phases. The classical BL equation gives a monotone 
water saturation profile for any Riemann problem; on the contrast, when the 
dispersive parameter is large enough, the MBL equation delivers non-monotone 
water saturation profile for certain Riemann problems as suggested by the ex- 
perimental observations. In this paper, we first show that the solution of the 
finite interval [0, L] boundary value problem converges to that of the half-line 
[0,-|-oo) boundary value problem for the MBL equation as L — -|-oo. This 
result provides a justification for the use of the finite interval boundary value 
problem in numerical studies for the half line problem. Furthermore, we extend 
the classical central schemes for the hyperbolic conservation laws to solve the 
MBL equation which is of pseudo-parabolic type. Numerical results confirm 
the existence of non-monotone water saturation profiles consisting of constant 
states separated by shocks. 
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4 1. Introduction 

5 The classical Buckley-Leverett (BL) equation [3 is a simple model for two-phase 

6 fluid flow in a porous medium. One application is secondary recovery by water-drive 

7 in oil reservoir simulation. In one space dimension the equation has the standard 

8 conservation form 



ut + {f{u))x = in 

(1.1) u{x,0)^0 

u(0, t) = ub 

9 with the flux function f{u) being defined as 

r u < 0, 

(1-2) 0<^<1 

I 1 u > 1. 



g = {(x,t) : x > 0,t > 0} 

X e (0, oo) 
t e [o,oo) 



10 In this content, m : Q — > [0, 1] denotes the water saturation (e.g. m = 1 means pure 

11 water, and u = means pure oil), its is a constant which indicates water saturation 

12 at a; = 0, and M > is the water/oil viscosity ratio. The classical BL equation 

13 (1.1 ) is a prototype for conservation laws with convex-concave flux functions. The 

14 graph of f{u) and f {u) with M = 2 is given in Figure [L1] 
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Figure 1.1. f{u) and f'{u) with M = 2. 



Due to the possibihty of the existence of shocks in the solution of the hyperbohc 



conservation laws (1.1), the weak solutions are sought. The function u G L°°(Q) is 



called a weak solution of the conservation laws (1.1) if 



Q 



M^+/(w)^^=0 foraU q^eC^iQ). 



Notice that the weak solution is not unique. Among the weak solutions, the entropy 
solution is physically relevant and unique. The weak solution that satisfies Oleinik 
entropy condition |19j 

(1.3) > s > for all u between ui and 

u — Ul u ^ Ur 

is the entropy solution, where u;, Ur are the function values to the left and right 
of the shock respectively, and the shock speed s satisfies Rankine-Hugoniot jump 
condition [TT l fTU ] 

(1.4) ^^/W-/K)^ 

Ul — Uj, 



The classical BL equation (1.1 ) with flux function f(u) as given in ( |1.2[ ) has been 
well studied (see [13] for an introduction). Let a be the solution of f'{u) ~ "^^i 
i.e., 

(1-5) « = 

^ ' V Af + 1 

The entropy solution of the classical BL equation can be classified into two cate- 
gories: 

(1) If < ub < oe, the entropy solution has a single shock at | = ^^^^^ ■ 

(2) If a < Us < 1, the entropy solution contains a rarefaction between ub and 



a for /'(ws) < I < f'{a) and a shock at 



/(") 



t 



These two types of solutions are shown in Figur e |1.2| for M = 2. In either case, 
the entropy solution of the classical BL equation ( |1.1[ ) is a non-increasing function 
of X at any given time t > 0. However, the experiments of two-phase flow in 
porous medium reveal complex infiltration profiles, which may involve overshoot, 
i.e., profiles may not be monotone [T. This suggests the need of modification to 



the classical BL equation (1.1 1. 
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Figure 1.2. The entropy so lutio n of the classical BL equation 



2, a 



0.8165). 



(a) 



consists of one shock at | = 



< uB ~ 0.7 < a, the solution 
a < uB = 0.98 < 1, the 



(b) 



solution consists of a rarefaction between ub and a for /'{ub) < 
f < /'(a) and a shock at f = 



(1.6) 







To better describe the infiltration profiles, we go back to the origins of (1.1 1. 
Let Si be the saturation of water/oil (i = w,o) and assume that the medium is 
completely saturated, i.e. S*^ + 6*0 = 1. The conservation of mass gives 

dt dx 

where ((> is the porosity of the medium (relative volume occupied by the pores) and 
Qi denotes the discharge of water/oil with + Qo = q, which is assumed to be a 
constant in space due to the complete saturation assumption. Throughout of this 
work, we consider it constant in time as well. By Darcy's law 



(1.7) 



ft 



Mi 



dx 



w, o 



where k denotes the absolute permeability, kri is the relative permeability and fii 
is the viscosity of water/oil. Instead of considering constant capillary pressure as 
adopted by the classical BL equation (1.1 1, Hassanizadeh and Gray [HI IS] have 
defined the dynamic capillary pressure as 

(1.8) Pc = Po - = PciSu,) - ' " 



dt 



where Pc{Sw) is the static capillary pressure and r is a positive constant, and 



dt 



is the dynamic effects. Using Corey [51[2n] expressio ns wi t h exp onent 2, kr^{Sw) — 



Q , rescaling x ^ 



for the water saturation u 
(1.9) 



X and combining (1.6)-(1.8), the single equation 



Sin is 



du 



d_ 
dx 



+ M(l - uf 



u- 



d_ 
dx 



02 k{l-ufu^ d fpc{u) 



(1 — u)2 + fJ-oU^ dx \ (j) 



du 
'~dt 



where M = ^ [35]. Linearizing the right hand side of (1.91 and rescaling the 
equation as in [3T1 [20] , the modified Buckley-Leverett equation (MBL) is derived 
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as 



(1.10) 



du 



dfju) 
dx 



eV 



dx'^dt 



where the water fractional flow function f{u) is given as in (1.2). Notice that, if Pc 
in (1.8) is taken to be constant, then (|1.9|) gives the classical BL equation; while 



if the dispersive parameter r is taken to be zero, then (1.10) gives the viscous BL 
equation, which still displays monotone water saturation profile. Thus, in addi- 
tion to the classical second order viscous term euxx, the MBL equation (1.10) is 



an extension involving a third order mixed derivative term e ruxxt- Van Dujin et 
al. showed that the value r is critical in determining the type of the solution 



profile. In particular, for certain Riemann problems, the solution profile of (1.10) 
is not monotone when r is larger than the threshold value , where was numer- 
ically determined to be 0.61 [21]. The non-monotonicity of the solution profile is 
consistent with the experimental observations ?]. 

The classical BL equation (1.1) is hyperbolic, and the numerical schemes for 
hyperbolic equations have been well developed (e.g. [TH [TSl |31 [SI UHl H] )■ The 



MBL equation ( 1.10 ), however, is pseudo-parabolic, we will illustrate how to extend 
the central schemes 181 [121 [13] to solve (1.10) numerically. Unlike the finite domain 



of dependence for the classical BL equation (1.1), the domain of dependence for 



the MBL equation (1.10) is infinite. This naturally raises the question for the 
choice of computational domain. To answer this question, we will first study the 
MBL equation equipped with two types of domains and corresponding boundary 
conditions. One is the half line boundary value problem 

ut + if{u))x = £Uxx + (?TUxxt in Q = {{x, t) : x > Q,t > 0} 
m(x, 0) — uq{x) X £ [0, oo) 

u(0,t) — gu{t), lim u{x,t) =0 te [0,oo) 

uq{0) = .gu(0) compatibility condition 

and the other one is finite interval boundary value problem 



(1.11) 



(1.12) 



+ {fiv))x = f^Vxx + (^"^TVxxt 

v{x,0) = vo{x) 
v{0,t) = g,(t), v{L,t) ^ h{t) 
= 5.(0), vq{L) = /i(0) 



Considering 
(1.13) 

Uq{x) = 



vq{x) for X G [0, L] 
for a; e [i, +oo) 



Q = {{x,t) : X e (0,L),t > 0} 

X £ [0, L] 

t e [0,oo) 

compatibility condition. 



gu{t)^gv{t) = 9{t), h{t) = 0, 



we will show the relation between the solutions of problems (1.11) and (1.12). To 



the best knowledge of the authors, there is no such study for MBL equation ( 1.10 ). 
Similar questions were answered for BBM equation [TJ |2] . 

The organization of this paper is as follows. Section [2] will bring forward the 
exact theory comparing the solutions of ( |1.11[ ) and ( 1.12| . The difference between 
the solutions of these two types of problems decays exponentially with respect to 
the length of the interval L for practically interesting initial profiles. This provides 
a theoretical justification for the choice of the computational domain. In section 
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[3) high order central schemes will be developed for MBL equation in finite interval 
domain. We provide a detailed derivation on how to extend the central schemes 



[TSIIT^ for conservation laws to solve the MBL equation ( 1.10 ). The idea of adopting 
numerical schemes originally designed for hyperbolic equations to pseudo-parabolic 
equations is not restricted to central type schemes only ([23l[24]). The numerical 
results in section |4] show that the water saturation profile strongly depends on the 
dispersive parameter r value as studied in [21\ . For t > t*, the MBL equation 
(1.10) gives non- monotone water saturation profiles for certain Riemann problems 
as suggested by experimental observations fZl. Section [s] gives the conclusion of the 
paper and the possible future directions. 



2. The half line problem versus the finite interval problem 



94 Let u{x, t) be the solution to the half line problem ( |1.11[ ), and let v{x, t) be the 

95 solution to the finite interval problem \1.12 1 . We consider the natural assumptions 

96 (1.131. The goal of this section is to develop an estimate of the difference between 



97 u and V on the spatial interval [0, L] at a given finite time t. The main result of 

98 this section is 







X e [0, Lo] 

X > Lq 



99 Theorem 2.1 (The main Theorem). If uq{x) satisfies 
(2.1) uo(x) = 

where Lq < L and Cu, are positive constants, then 

II u{-,t) - v{-,t) 11^1 < Di,,^r{t)er^ +D2..eAty 

for some < A < 1, Di-^^rit) > and -D2:e,r(^) > 0, where 



iin.oiim 



Y{x,tY + {e^Y^{x,t)Ydx 



108 
109 
110 



112 
113 
114 



Notice that the initial condition (2.1 1 we considered is the Riemann problem. 
Theorem 2.1 shows that the solution to the half line problem ( |1.11| can be approx- 
imated as accurately as one wants by the solution to the finite interval problem 
(1.12) in the sense that Z^i.^ ,-(t), D2-e,Tit), and ^'^^^7^°^ can be controlled. 



105 To prove theorem 2.1 we first derive the implicit solution formulae for the half 



line problem and the finite interval problem in section 2.1 and section 2.2 respec- 
tively. The implicit solution formulae are in integral form, which are derived by 
separating the x-derivative from the t-derivative, and formally solving a first order 
linear ODE in t and a second order non-homogeneous ODE in x. In section |2.3[ 
we use Gronwall's inequality multiple times to obtain the desired result in theorem 

nm 

2.1. Half line problem. In this section, we derive the implicit solution formula 
for the half line problem ( |1.11[ ) (with gu(t) — g{t)). To solve ( |1.11| , we first rewrite 
(1.11) by separating the .T-derivative from the t-derivative. 



(2.2) 



dx'^ 



Ut H u 



er 
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115 By using integrating factor method, we formally integrate (2.2) over [0,t] to obtain 



(2.3) 



117 Furthermore, we let 
(2.4) 



u — e " Mo = 



1 



ifiu))x e 



ds. 



A 



118 then (2.31 can be written as 

ft 



(2.5) A" - 



1 



A = 







ds, 



where ' — 



dx ' 



119 Notice that (2.5) is a second-order non-homogeneous ODE in s-variable along with 

120 the boundary conditions 

A{Q,t)=u{Q,t)-e-^ua{Q) = g{t) ^ ^ g{Q) , 
A(oo,t) = u{oo,t) — e^"" uq{go) = 0. 



(2.6) 



121 To solve ( |2.5[ ), we first solve the corresponding linear homogeneous equation with 

122 the non-zero boundary conditions (2.6). We then find a particular solution for the 

123 non-homogeneous equation with zero boundary conditions by introducing a Green's 

124 function G{x, f ) and a kernel K{x,C) for the non-homogeneous terms u and {f{u))x 

125 respectively. Combining the solutions for the two non-homogeneous terms and the 

126 homogeneous part with boundary conditions, we get the solution for equation (2.5) 

127 satisfying the boundary conditions (2.6): 



A{x,t) 



1 



(2.7) 



1 



^0 



G{x,0<i,s)t 



ds 



-\-oo 



K{x,i)f{u)e-— d^ds 
[g{t)-e-^g{0)^e-^ 



128 where the Green's function G{x,^) and the kernel K(x,S^) are 
(2.8) G{x,0 = ^(e-^-e- — 



(2.9) 



K{x,0 



dG{x,0 



I I e 



sgn(a::-Oe' 



129 To recover the solution for the half line problem ( |1.11[ ) , we refer to the definition of 

130 A in (2.4). Thus, the implicit solution formula for the half line problem ( |1.11[ ) is 



u(x, t) 



+ 00 



e '-^ — e 



m(^, s)e " dS^ds 



(2.10) 



+ 



2eV Jo Jo 



+00 



e + sgn{x — S^)e \ f(u)e 



ff(0) ) e 



dS, ds 



'uo{x). 
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131 2.2. Finite interval problem. The implicit solution for the finite interval problem 

132 ( 1.12 1 (with gv{t) = g{t)) can be solved in a similar way. The only diflFerence is that 

133 the additional boundary condition h(t) at x = i in (1.12) gives different boundary 

134 conditions for the non-homogeneous ODE in x-variable. Denote 

(2.11) ^v-e^^vo, 

135 then it satisfies 



L\// 



(2.12) (A^) 



1 



-A' 



1 



1 



{f{'v))x I e " cis where ' 



d_ 
dx 



with the boundary conditions 

A'^{Q,t)^v{{),t)-e-^VQ{Q)=g{t)-e-^g{Q), 
A^{L,t) = v{L,t) -e'^voiL) ^ h{t) - e- ^ h{0) . 
136 These boundary conditions affect both the homogeneous solution and the par- 



137 ticular solution of (2.12) as follows 

1 

1 



A^{x,t) = - 



(2.13) 



Jo 

* r K\x,Of{v)e-'^d^ds 
Jq 

Cl{t)(f>l{x) + C2{t)(j)2{x) 



138 where the Green's function G^{x,£,), the kernel K^{x,£,) and the bases for the 

139 homogeneous solutions are 



(2.14) 

(2.15) 

(2.16) 
(2.17) 



= If^ (e^ 

2(e^ - 1) ^ 



e — e — e 



K\x,£,) = - 



-K 2I.-(x + {) 

gjv^ — e 



2(e'v^ - 1) 



-|-sgn(x — ^)e — sgn(a; — f)e 



'/'lis:) = ^ 



C2(i) = /t(t)-e-"/i(0), 



and 02(2;) 



141 Thus, the implicit solution formula for the finite interval problem (1.12) is 

v{x, t) = — 



2t^T^{e'^ - 1) Jo Jo 



f x+{ 2i,-(x+i;) Ix-e 



(2.18) 



2g2^(g,^ - 1) Jo Jo 



+ sgn(x — ^)e 



-sgn(a; — ^)e '-^ ) /(w)e " ds 



ci(i)0i(a;) + C2(i)02(a;) + e "Wo(a;) 
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2.3. Comparisons. In this section, we will prove that the solution u{x, t) to the 
half line problem can be approximated as accurately as one wants by the solution 



v{x,t) to the finite interval problem as stated in Theorem 2.1 

Due to the difference in the integration domains, we do not use ( |2.10 1 and (2.18) 



146 directly for the comparison. Instead, we decompose u{x,t) {v{x^t) respectively) 

147 into two parts: U{x,t) and UL{x,t) {V{x,t) and VL{x,t) respectively), such that 

148 U{x,t) {V{x,t) respectively) enjoys zero initial condition and boundary conditions 

149 at a; = and x = L. We estimate the difference between u{-,t) and v{-,t) by 

150 estimating the differences between UL{-,t) and fL(-,t), U{-,t) and V{-,t), then 

151 applying the triangle inequality. 



2.3.1. Definitions and lemmas. To assist the proof of Theorem |2.1| in section ["2.3.3[ 
we introduce some new notations in this section. We first decompose u{x, t) as sum 
of two terms U{x,t) and UL{x,t), such that 

u{x,t) ^ U{x,t) + UL{x,t) x e [0, +oo) 

153 where 

(2.19) UL = er^uo{x) + ci{t)e~ 



+ {u{L,t) ~ cx{t)e 



e " uq{V)\ <i)i{x) 



154 and c\(i) and (^^{x) are given in (2.16) and (2.17) respectively. With this definition, 

155 takes care of the initial condition uq{x) and boundary conditions g(£) at x = 

156 and X = L for u(x,t). Then U satisfies an equation slightly different from the 

157 equation u satisfies in ( [l.lip : 

(2.20) 

Ut - eUxx - e^rUxxt = {ut - eu^x - e^ruxxt) - (("l)* - (-[ul)xx - e^T{uL)xxt) 
= -{f{u))^ + -UL{x,t) 

158 In addition, U{x,t) has zero initial condition and boundary conditions at x = 

159 and X ^ L, i.e., 

(2.21) U{x,0)^0, C/(0,t) = 0, U{L,t)=0. 
Similarly, for v{x,t), let 

v{x,t) ^V{x,t) +VL{x,t) xe[0,L] 

160 where 

(2.22) VL = e-^voix) + ci(t)0i(x) + C2(t)</)2(a;) 

161 and ci(t), C2{t) and (j)i{x), 4>2{x) are given in (2.16) and ( |2.17 ) respectively. With 

162 this definition, takes care of the initial condition vq{x) and boundary conditions 

163 g{t) and h{t) at a; = and x = L for v{x,t). Then V satisfies an equation slightly 

164 different from the equation v satisfies in ( |1.12 ): 



(2.23) 
165 with 



Vt - eVx._ 



e^TVxxt^~{f{v)),, + —VL{x,t) 



(2.24) 



V{x, 0) = 0, V^(0, t) = 0, V{L, t) = 0. 
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Since, in the end, we want to study the difference between U{x,t) and V{x,t), we 
define 

W{x,t) = V{x,t)-U{x,t) for xe[0,L]. 

166 Because of (2.20) and ( |2.23 1, we have 

(2.25) Wt - eW^^ - e^TW^^t - - (/(«) - + —{vl - ul). 

er 

167 In lieu of (2.211 and ( 2.24[ ), W{x,t) also has zero initial condition and boundary 

168 conditions at a; = and x ~ i.e., 

(2.26) W^(a;,0) = 0, W{Q,t) = Q, W{L,t)=0. 

169 Now, to estimate — we can estimate — \\V — U \\ and estimate 



170 II Ul — Vl II separately. These estimates are done in section 2.3.3 



171 Next, we state the lemmas needed in the proof of Theorem |2.1[ The proof of the 

172 lemmas can be found in the appendix [A| and [22 . In all the lemmas, we assume 

173 < A < 1 and Uf){x) satisfies 



(2.27) 



uq{x) 



Cu x£ [0, io] 

x> Lq 



174 where Lq < L and C„ are positive constants. Notice that the constraint A G (0, 1) 

175 is crucial in Lemmas 12.31 12.41 



176 Lemma 2.2. f{u) = 



2 + M(l— u)5 



< Du where D 



fM 



a 



and a 



M 
M+1 ■ 



177 Lemma 2.3. (i) 

178 (a) J^°° 



+ 00 



5±5 



e "v^ — e 



e — e '-^ 



e — e 



A2 ■ 



" < e(l-A) ■ 



Ain 



180 Lemma 2.4. ('ij /q^°° e + sgn(x — ^)e '-^ e < ■ 



(in) J^°° 



e + sgn(a: — ^)e ^■■^ 
e + sgn(a; — ^)e 



A:r-g 



e(l-A) 



e=v^|uo(0|de < 2C„eV^e^ 



183 Lemma 2.5. ('ij 0i(a:) — e "■■^ —e \4>2ix)\ . 

184 (ii) \(j)2{x)\ <lforxe [0,i] . 

185 (lii) \(l3'2{x)\ < ^ j/e < 1 /or X e [0, L] . 

186 Last but not least, the norm that we will use in Theorem |2.1| and its proof is 



(2.28) 



Y{x,t)^ + {eV^Y^{x,t)ydx. 



187 2.3.2. A proposition. In this section, we will give a critical estimate, which is es- 

188 sential in the calculation of maximum difference || Ui(-,t) — VL{-,t) ||^ in section 

189 2.3.3 By comparing UL{x,t) and VL{x,t) given in (2.191 and (2.22) respectively. 



190 it is clear that the coefficient u{L,t) — Ci(t)e '-^ — e '^Uq{L) for (f)2{x) appeared 
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191 in (2.19) needs to be compared with the corresponding coefficient C2{t) for 4>2{x) 

192 appeared in (2.22). We thus define a space-dependent function 

(2.29) Ucjx^t) = u{x,t) - ci{t)e~^ - e'^uoix) 

193 and estabfish the foUowing proposition 
Proposition 2.6. 

b-rt AL t (b-r-l)t A(I,-I,o) 

(2.30) \Uc2{L,t)\ < ar{t)e~e +Cr — e^^e 

194 for some parameter- dependent constants a^, and Ct- 



195 Proof. Based on the impficit solution formula (2.10) derived in section 2.1 Lemma 

196 2.2 and the relationship between 17^2 and u given in (2.29), we can get an inequality 



197 in terms of Uc. 



\UcJx,t)\ < 



1 



t n + OO 



t f + OO 



JO 



e — e 



\U,,{ts)\e-^ d^ds 



"'0 



e 



\ci{s)\e~ 



d^ ds 



(2.31) 



D 

2e^ 



e =^ — e 

t /'+00 



JQ 



JO 

t /' + 00 



e + sgn(a; — ^)e \Uc2iC,s)\e " d^ds 



e + sgn(a; — ^)e |ci(s)|e 'v^e " ds 



e + sgn(x — ^)e |uo(^)|e d^ds 



lo Jo 

198 To show that Uc2{x,t) decays exponentially with respect to x, we pull out an 

199 exponential term by writing Uc2ix,t) = e e~^U(x,t), where < A < 1, such 

200 that 



(2.32) 



U{x,t) = e-^-^e"Uc2{x,t), 



201 then (2.31 ) can be rewritten in terms of U{x, t) as follows 



ty{x,t) 



< 



1 



t /' + 00 



t p+oo 
^0 



e «^ — e 



ds 



Jo 

t /'4-00 



|ci(s)| e 'v^ ds 



(2.33) 



Jo 



e-^ \uoiO\ d^ds 



D 

2e^ 



t p-\-oo 



^0 

t n-\-oo 



U(ts] 



d^ ds 



Jo 



Jo 



_£±5 la:-!; Ax-{ ^ 

e + sgn(x — ^)e |ci(s)| e d^ds 



e + sgn(a; — ^)e e |uo(^)| d^ds 
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202 Because of Lemmas 

203 on ( |2.33[ ) : 



2.3 



2.4 



we can get the following estimate for 



Ui;t) 



based 



< 



1 - A2 7o 



\U{-,s)\oods + 



e(l - A) Jo 



+2C^,e^/7e«^/^ / Ids 



(2.34) 



D 
2^ 



2e^ 
1 - 7o 



\Ui-,s)\^ds + ey/^ [ 1 + 



e(i-A);yo 



|ci(s)|e«T ds 



\ci{s)\e'^ ds 



+2C^,e^/Te-v^ / Ids 



< / —\U{-,s)\oods + 
er 







ds 







er 



where 



br 



1-A2 ' 
|ci(-)U(l + ^\/^(e(l- A) + l)) 



2e(l - A) 



By Gronwall's inequality, inequality (|2.34|) gives that 
U{-,t) 



< 



oo 



dr{t — s) b^jt-s) 







er 



e " ds < [ fli-e'^ + Cr — e^-^ e " 



er 



204 Hence < V{-,t) 

205 Uc.2{x,t) decays exponentially with respect to x. In particular, when x = L, we 



oo 



e«v^e »T < I a^gET + c,- — e I e =^ e e " i.e., 



206 have 
(2.35) 



|(yc9(-t^,i)| < Ore " e -\- Cr — e " e 

er 



207 as given in ( 2.30 1. 



□ 



208 2.3.3. Proof of Theorem \2.1\ In this section, we will first find the maximum dif- 

209 ference of || i) — i) || , then we will derive || ul(-, t) — VL{-,t) and 

L .E ,T 

210 II W{-,t) 11^1 ~ II U{-,t) — V{-,t) 11^1 . Combining these two, we will get an 

211 estimate for II u(-, t) — !;(•, t) II rri 

Proposition 2.7. If ua{x) satisfies 12.21 ), then 

II ul-vl\\^< Ei,,^r{t)e'^ + E2.,,^r{t)e'^^^^ 

212 where Ei-^t^T-it) = \ci{-)\oo + o-T^^ and i?2;£,T(i) — Cr^e' ' . 

Proof. By the definition of and given in ( |2.19 1 and (2.22) and the assumption 
that uo(a^) = ^0(2;) for x € [0, L], we can get their difference 



UL{,x,t)-VL{x,t) = ci{t) [e -(t>i{x)j + (^C/c,(L,t) -/i(t) + e-^^/i(0)j <?!)2(a;) 
213 Combining Lemmas 2.5[ il), [275| pll), inequality (2.35), and h{t) = 0, we have 

(2.36) \\uL{-,t)-VL{-,t)\\^ < Si.,,,(Oe~^ +^2;e,r(0e"^^^ 
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214 where 



(2.37) Ei., ^{t) = |ci(-)|oo + a^e and ^^2;£,r(i) = c^ — e* ' . 



□ 



Proposition 2.8. If uq{x) satisfies {2.21), and Ei-^^rit), E2-e,Tit) are as in propo- 
sition\2. 71 then 



XL \(L-La) 

\\uL{-,t)~VL{-,t)\\Hi < VbL [ Ei,,^rit)e + E2;e,r{i)e 



217 Proof. Because of the definition of ul and given in (2.191 and (2.22), Lemma 

218 2.5 iii| and inequaUty (2.35), we have that 



(2.38) 



iuLi;t) - VLi;t)). \\^ < \ciit)\e '-^ |0^(a:)| + \U,,{L,t)\ |</>^(a:)| 

2 / , , ^AL_ , , KL-Lq) 



219 Now, combining (2.36) and (2.38), we obtain that 



(2.39) 



\\uL{-,t)~VL{-,t)\\^i 



<V5L (Ei.,,^r{t)e + S2;e,r(t)e 



□ 



Proposition 2.9. If uq{x) satisfies ^2.27 ), then 



XL X(L-Lq) 

||W^(-,i)|lm <7l;..r(t)e +72;.,r(t)e ^v^^ 



221 where the coefficients are given by 



7l;e,r(0 = 6 ^«'v^ 

(M + 1)^1 

(2.40) 72;e,r(0 = e 2Af<v^ 



V 2Af 

(i£±ll^ /(M + l)Vr 

V 2M 

t (bT-l)t 

e " 



l)^/L(^|cl(.)U + J(e^-l) 



eT(6, - 1) 



1 j 

1 

(67^ 



(e " - 1) 



223 Proof. Multiplying the governing equation of W (2.25) by 214^, integrating over 

224 [0, i], and using integration by parts, we get 



/ W'^ + {t^W^fdx 
Jo 



d^ '■^ 
dt 



2W^ dx + I 2W^ ifiv) - f{u)) dx+ — I W{vL - ul) dx. 
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225 Therefore, using the norm we defined earlier in (2.281, and f'{u) < 

226 we have 



(M+l)^ _ ^ 



< 2 I \W,\\f'(,^)\\v-u\dx+^\\vL-UL\\^\\Wi;t)\\j,,^^^_^ 

< 2C I m\i\W\ + \\vL~UL\Udx+^\\vL-UL\\^\\Wi;t)\\Hr 







< 



2C 



(\\W{;t)\\m +\\vL~UL\LVL\\W(;t)\\^^ ) 



2C 



II Wi;t) \\' + ("^ + 1 W II «i - IL II Wi;t) II . . 



Hence 
d 



C 



c i_ 



\/i II - 



227 By Gronwall's inequality and (2.36) 



< 



t 



Ct 1 


' C 






< e"^ 1 




Vl f 











ds 



< 



( C \ 



< e 
+e 
Hence 



/ C f 

+ _ 

V e^/T er 

/ C f 

+ _ 

ct / C 1 



i|ci(-)| 



(e~ - 1) e 'v^ 



e-v/r er / er V 6,- — 1 



for 



" -l))e 



A(I--I.o) 



Iim-,0llm <7l;.,r(t)e +72;e,r(i)e" 



228 where 7i;e,T(0 ^-^id "f2:e.T{t) are given in (2.40). 

229 Now we are in the position to prove the main theorem of this section. 

230 Theorem 2.10. If uo{x) satisfies 

uo{x) -- 



□ 







X e [0, Lo] 

X > Lq 
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231 where Lq < L andCu, are positive constants, and Ei-^^T{t) , E2-e,T{t) , Ji-a^rit) , "f2;e,Tit) 

232 are as in (2.371 and (2.401 , then 



(2.41) \\u{;t)-v{-,t)\\^,^^^<Di,,,r{t)e '-^ + D2-e^r{t)e 

for some < A < 1, and 

Due.rit) = 7l;e,r(t) + V5LEi,,^r{t) , -D2;e,rW = 72;e,r(t) + V5LE2;e,r{t) ■ 

Proof of the Main Theorem. 



II ^i-^t) - v{;t) WhI , _ < II W{;t) 11^; _ + II VL{;t) - ULi-,t) ||^i^ 



' L.C.T 

XL 



where 



V 2M 

+ \/5Z(|c(-)|oo + flrC^), 
D2;eAt) ^l2-cAt) + V5LE2,e.r{t) 



I ^ ^ + 1 ) x/I ( ^|ci(-)|oo + ^(e^ - 1) 



hr 



t (bx-ljt 1 (bT-l)t 
P ex — (p e-r — \ 



r— t (fcx-l)* 

'oLcr — e " 
er 



□ 



234 This rcsuh gives that || t) — || exponentially delays in L. This 

235 theorem shows that if and ^'^^^^"'^ converge to infinity, then the solution 

236 v{x,t) of the finite interval problem converges to the solution u{x,t) of the half 

237 line problem in the sense of || • ||^i . This can be achieved either by letting 

238 L — (X) or e — > 0. For example, in the extreme case, e = 0, the half line problem 



239 ( 1.11 1 becomes hyperbolic and the domain of dependence is finite, so, certainly, one 



240 only need to consider the finite interval problem. This is consistent with the main 

241 theorem in the sense that for a fixed final time t, if XL > brt and A(L — Lq) > 

242 {br-l)t, i.e., L > max( ^^^r^), then || u{-,t) - v{-,t) ||^i < Di,^^r{t)e^^ + 

X(L-La) 

243 £'2;e,r(0e 



2.10 



gives a theoretical justification for 



as e — > 0. Theorem 

244 using the solution of the finite interval problem to approximate the solution of the 

245 half line problem with appropriate choice of L and e. Hence in the next chapter. 



246 the numerical scheme designed to solve the MBL equation (1.10) is given for finite 

247 interval problem. 
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3. Numerical schemes 



249 To numerically solve the MBL equation ( 1.10 1, We first collect all the terms with 



258 
259 
260 
261 
262 
263 
264 



250 time derivative and rewrite MBL equation (1.10) as 
(3.1) 



c)t + {f{u))x = 



251 By letting 
(3.2) 



W = U ~ e TU.r 



252 MBL equation (3.11 can be written as 



(3.3) 



Wt + {f{u))x = 



253 Now, the new form of MBL equation (3.3) can be viewed as a PDE in terms of w, 



254 and the occurrence of u can be recovered by (3.2). Equation (3.3) can be formally 

255 viewed as 



(3.4) 



wt + ifiil- 



256 which is a balance law in term of w. We adopt numerical schemes originally designed 



257 for hyperbolic equations to solve the MBL equation (3.1), which is of pseudo- 



parabolic type. The local discontinuous Galerkin method has been applied to solve 
equations involving mixed derivatives Uxxt term j23l I24j . To the best knowledge of 
the authors, the central schemes have not been applied to solve equations of this 
kind. The main advantage of the central schemes is the simplicity, "the direction of 
the wind" is not required to be identified, and hence the field- by-field decomposition 
can be avoided. In this chapter, we demonstrate how to apply the central schemes 



to solve the MBL equation (3.1 ) 



265 
266 
267 



3.1. Second-order schemes. In this section, we show how to apply the classical 
second order central schemes [T8j originally designed for hyperbolic conservation 



laws to numerically solve the MBL equation (1.10 1, which is of pseudo-parabolic 
type. To solve (3.3), we modify the central scheme given in jl8j. As in [T5], at each 
time level, we first reconstruct a piecewise linear approximation of the form 



(3.5) L,{x,t) = w,{t)^{x-Xj)^, 



Xj^l <X<X^+l. 



270 Second-order accuracy is guaranteed if the so-called vector of numerical derivative 

271 which will be given later, satisfies 



(3.6) 



w'^ dw{xj,t) 



3 



dx 



0{Ax) 



272 We denote the staggered piecewise-constant functions Wjj^i(t) as 



(3.7) 



Ax 



w{x, t) dx. 
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273 Evolve the piecewise linear interplant (3.5) by integrating (3.3 1 over [xj,Xj+i] x 

274 [t,t + At] 



(3.8) 



1 
Ax 

e 

Ax 



t+At 



f{u{xj+i,s))ds 

*+'^* d^u{x, s) 

dx'^ 



t+At 



f{u{xj,s))ds 



dx ds 



275 We calculate each term on the right hand side of (3.8 ) below. For w^^i (t), applying 

276 the definition of Lj{x,t) and given in (3.5) to (3.7), we have that 



(3.9) 



Ax 



Lj{x, t) dx + 



Ax 



Lj+i{x, t) dx 



3 3- 

277 The middle two integrands can be approximated by the midpoint rule 



(3.10) 



t+At A , 

/(u(x„ s)) ds = f{u{x,,t + — ))At + 0{At') 



t+At 



At, 



f{u{xj+i,s)) ds = fiu{xj+i,t + ^))Ai + O(At^) 



if the CFL condition 



A • max 



df{u{w{x,t))) 
dw 



< 



1 



where A 



At 
Ax 



is met. For MBL equation (3.3), we have that at < > 0, 

u — €^TUxx=w, u{0)=w{0), u{L)~w{L). 



Let v{x) = iL-^)^io)+-^W ^ then 



i{x) = [Iw]{x) ^ v{x) + - / [wiy) -v{y)]K{x,y)dy 



where 



fc=l 



Hence the eigenvalues for / are 



A*. 



1 



l + (^)2e2r 



< 1, 



k = 1.2,3. 



Therefore, the CFL condition is 



At 

— — • max 

Ax Xj<x<Xj^i 



dfiu{w{x,t))) 



dw 



At 

— — • max 

Ax Xj<x<Xj^i 
fe=l,2,3... 



dfiuix,t)) 



du 



• Afc < ^ • 2.2 < 
Ax 
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278 In the numerical computations in chapter |4j we chose ^ — 0.1. In (3.10), to 

279 estimate + ^)'s, we use Taylor expansion and the conservation law (3.3 1: 

At, , . dw At ^ o. 



(3.11) 



, , , d'^u df .At , , n, 
w,{t) + {eAxD^u,-f')^, 



where D is the discrete central difference operator 



— 2uj + Uj^i 



280 and the second-order accuracy is met if 

fj _ df{u{xj,t)) 



(3.12) 



Aa; 



dx 



0{Ax). 



281 The choices for {w^} in (3.6) and {/j} in (3.12) can be found in [T8', and we chose 
(3.13) w'^ = AfAf { Au;^.+ 1 , Aw^_.} , /j - Af M{A/^.+ 1 , A/^_ i } 

282 where MM{x,y} = minmod(a;,y) = i(sgn(a;)+sgn(y)) -Minda;!, |y|) and Awj_^i = 

283 Wj^i —Wj. Combining (3.8)-(3.10), we obtain 



w^+i it + At) ^w^+ijt) 



(3.14) 



-X[f{u,^,it+^)^f{u,{t+^))] 



e 
Ax 



dx"^ 



dx ds 



Next, we will re-write (3.14) in terms of u. {uxj.),_^i is approximated as 



{^xx ) j 



1 

Ax 



Ax 



{Ux{Xj + l,t) - U^{Xj,t)), 



284 and using the cell averages, it becomes 



(3.15) 



Ax V Aa; Ax 

_ Uj+3/2 - 2Mj+l/2 + •»j-l/2 

^ (Aa;)2 

2 



Notice that the linear interpolation (similar to (3.5|) 



Lj^i {x,t + At) = u^^i {t + At) + (x - xj^i) 



5' Aa 



for Xj < X < Xj^i 



and the cell average definition (similar to (3.7)) 



u^+iit + At) 



At 



^3 + 1 



u(x, t + At) dx 



ensure that 



w,+ i(< + Ai) -Wj+i(t + At), 
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285 and the convertion between u and w is done using the following relation 
(3.16) {I-e^TD^)u^w. 



286 Hence re-writting (3.14) in terms of u gives the staggered central scheme 



(3.17) 



{I- 






At 


A[/(u,+i(i+ — 


e 


pt-{-At pXj 


Ax 


J t J Xj 



2 " 
dx ds 



We will focus on the last integral in (3.17). There are many ways to numerically 
calculate this integral. We will show two ways to do this in the following two 
subsections, both of them achieve second order accuracy. 



290 3.1.1. Trapezoid Scheme. In this scheme, we use the notion (3.7) and the trapezoid 

291 rule to calculate the integral numerically as follows: 



(3.18) 



*+'^* d'^u{x,s) 
dx"^ 



t+At 



dx ds — Ax I {uxx)j+^{s) ds 



AxAt 



292 with 0{At^) error. Combining with (3.15) and (3.17), we can get the trapezoid 

293 scheme 



(3.19) 



I-{e'r+'-^)D-^]u^^.{t + At)=[l-{e^r- 



eAt, 



-A 



294 The flow chart of the trapezoid scheme is given in (3.20) 



(3.20) 



\3.16\ 

Uj [t) > Wj [t) 



u;,(i+f ) 



ii,+ i(t + Ai) 



3.1.2. Midpoint Scheme. In this scheme, we use the notion (3.7) and the midpoint 
rule to calculate the integral numerically as follows: 



^t+Zit p 
Jt Jx 



dx"^ 



dx ds — Aa 



t+At 



K^)i+i(s) ds 



At, 



= AxAt{u^)^+i{t+—) 

295 Combining with ( 3.15[ ) and (3.17), we can get the midpoint scheme 
{I-e\D^)u^+^{t + At) =z2),+ i(t) 



(3.21) 



A[/(%+i(t+^) -/(«,(*+ ^))] 
eA<i?2y (;+ ) 
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296 The flow chart of the midpoint scheme is given in (3.22) 
(3.22) 

i3S\ ' 13.21 




297 3.2. A third order semi-discrete scheme. Similarly, we can extend the third 



298 order scheme to solve MBL equation (1.10), however, it is more involved. But 



299 the third order semi-discrete central scheme proposed in |12| can be extended to 

300 solve the MBL equation in a straightforward manner. In order to make the paper 

301 self-contained, we include the formulation below. 



dt ~ Ax 



eQ,{t) 



302 where w{x,t) denotes the cell average of w 



Wj{t) = -— / w(x,t)dx, 



303 Hj^i/2{t) is the numerical convection flux and Qj{t) is a high-order approximation 

304 to the diffusion term u^x- 



./■(<+i/2(0) + /K+i/2(0) a,+i/2(t) 



'^'j>l/2W-^j + l/2W 



where 'uj_^_^^^{t),u'^^^^^{t) denote the left and right intermediate values of u{x,t") 
at Xj^i/2, and their values are converted from the 'wj_f_j^^2(^) ^ "^^+1/2^^) using (3.2). 
The way to calculate wj^^^^{t), 'w^_^_l/2{^) and aj^i/2{t) is 



%+i/2^*) - ^J+i 2" ^ 8 — ^1+^' 

%+i/2W = + -Y^i + 

( df _ df 
aj+i/2(<) = max I ^(^^7+i/2(^)), ^(4+i/2(*)) 
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321 
322 



324 
325 



329 
330 



where 



A, 

CL 

ISl 



1 

Ax 

2wc- 
ai 



wc 



2w^ 



wr{w' 



j+1 - w^)+wc- 
- 2wl 



■ w 



i+1 



1/4, 



IS, 



c 



y 

CR = 



13, _„ 



Ax2 

a, = 



Ci 



w 



cc = 1/2, eo = 



e{C,i?,L} 



10" 



n\2 



2*" 



305 The diffusion Uxx is approximated using the following fourth-order central differ- 

306 encing form 



(3.23) 



Q,{t) 



307 The unique feature of this scheme is that the discretization is done in space first, and 

308 then the time evolution equation can be solved as a system of ordinary differential 

309 equations using any ODE solver of third order or higher. In this paper, we simply 

310 use the standard fourth order Runge-Kutta methods. Notice that to achieve the 



311 third order accuracy, the linear solver that converts u from w using (3.2 1 need also 



312 to be high order, and (3.23) is used to discretize u^x in our convertion. 



313 4. Computational results 

314 In this section, we show the numerical solutions to the MBL equation 

(4.1) Ut + {j{u))x^ tUxx^ i^TUxxt 

315 with the initial condition 

ub if x = 



(4.2) 







if x > 



316 and the Dirichlet boundary condition. 

317 Numerically, it is not practical to solve the half line problem (4.2), and one 



has to choose an appropriate computational domain. Theorem 2.10 in Chapter [2] 
provides a theoretical bound for the difference between the solution to the half line 
problem and that to the finite interval problem. However, the estimate ( 2.41[ ) in 
Theorem 2.10 includes time-dependent parameters D^-^^ and D2-e,T{t), which 

Therefore, we numerically demonstrate how the 

We choose t = 5, ub = a = 



cannot be obtained analyticaly. 
computational domain size affects the solution 



and e = 0.001 as an example here. Figure 4.1 shows the snapshot of the solutions at 
t = 0.1, t = 0.5 and t = 1 for computational domain [0, L] with L = 0.25, L = 0.75 
and L = 1.25.^ 

t = 0.1, the leading shock is located at ^^^^ x 0.1 = 1.02 x 
- 0.25, L = 0.75, L — 1.25 all exceed the leading shock location. 



In Figure 
0.1 = 0.102, aEi 



Hence all the three computational domains deliver visually indistinguishable results. 
Whereas, in Figure 4.] |[b)[ t — 0.5, the leading shock is located at 1.02 x 0.5 = 0.51, 
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(a) t = 0.1 



(b) t = 0.5 



(c) i = 1 



1 






1 


0.8 






0.8 


0.6 






0.6 


0.4 






0.4 


0.2 






0.2 













0.1 



0.2 



0.2 0.4 0.6 




(a) t = 0.1,^ 



Figure 4.1. Numerical solutions of MBL (4.1) at 
t = 0.5, |(c)| ^ = 1 using the trapezoid scheme (3.20). ' — ' ' — 
' denote the numerical solutions corresponding to computational 
domain [0,L] with L — 0.25, L = 0.75 and L — 1.25 respectively. 

The parameter values are t — 5, 
Ax ^ ^, At = O.lAx. 



ub 



e = 0.001, 



331 L — 0.25 is shorter than the computational domain needed to capture this shock, 

332 hence the numerical solution halts at x — 0.25. On the contrast, L — 0.75 and L = 

333 1.25 are both large enough to capture this shock front. Similarly, in Figure 4.][c) 



334 t = 1, the leading shock is located at 1.02. L = 0.25 < 1.02 and L = 0.75 < 1.02 

335 both result in wrong solution profiles. More specifically, both solutions halt at the 

336 boundary of the insufficient computational domain. But L = 1.25 > 1.02 is large 

337 enough to capture the correct solution profile. 

In the rest of this chapter, all the computational domains [0, L] are therefore 
chosen based on the principle: 

L > leading shock speed x computational time. 

338 In addition, numerical solutions for larger L's, for example, L = 1.75, L = 2.5, 

339 L = 5, L = 10 are also sought. For all these larger L's, the numerical solutions 

340 are all consistent with that corresponding to L = 1.25 up to i = 1. This confirms 

341 that it is not necessary to take L too much larger than leading shock speed x 

342 computational time. 

343 To validate the order analysis given in chapter [3] for various schemes proposed, 

344 we first test the order of our schemes numerically with a smooth initial condition 

uq{x) = ubH{x — 5, 5), 

345 where 

r 1 if a; < 

H{^,0^{ l-i(l + f + ^sin(f )) if ^(<x<^ . 
[ if a; > ^ 

346 The final time T = 1 was employed, so that there was no shock created, e in the 



347 MBL equation (4.1 1 is taken to be 1, M is taken to be 2, and the computational 

348 interval is [—10,20]. The Li,L2,Loo order tests of the trapezoid scheme and the 

349 third order semi-discrete scheme with different parameter r value and the initial 

350 condition ub are given in Tables |4.1[ |4.2| Table |4.1| shows that the trapezoid rule 

351 achieved second order accuracy for all the tested cases in Li,L2,Lac, sense. Table 
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IN 


II 

II 2 


order 

1 


II 

^Ax — M Ax 
II 2 


order 

2 


II 

^Ax — M Ax 
II 2 


order 






60 


7.5416e-03 


- 


2.5388e-03 


- 


1.5960e-03 


- 




= 0.9 


120 


1.9684e-03 


1.9379 


6.7288C-04 


1.9157 


4.4066e-04 


1.8568 




0.2 


240 


4.9891e-04 


1.9802 


1.7645e-04 


1.9311 


1.2529e-04 


1.8144 






480 


1.2589C-04 


1.9865 


4.5366e-05 


1.9596 


3.3205e-05 


1.9158 






60 


8.0141C-03 


- 


2.6069e-03 


- 


1.4989e-03 


- 




= 0.9 


120 


2.1502C-03 


1.8981 


7.0452e-04 


1.8876 


4.2221e-04 


1.8279 




1 


240 


5.5697e-04 


1.9488 


1.8259e-04 


1.9480 


1.1283e-04 


1.9038 






480 


1.4104e-04 


1.9815 


4.6109e-05 


1.9855 


2.8719e-05 


1.9740 






60 


1.3102e-02 


- 


4.1784e-03 


- 


2.2411e-03 


- 




= 0.9 


120 


3.6201e-03 


1.8557 


1.0994e-03 


1.9263 


6.1060e-04 


1.8759 




5 


240 


9.6737e-04 


1.9039 


2.8089e-04 


1.9686 


1.5667e-04 


1.9625 






480 


2.5825e-04 


1.9053 


7.1250e-05 


1.9790 


3.9286e-05 


1.9956 






60 


6.4427e-03 


- 


2.1578e-03 


- 


1.1682e-03 


- 




— fj 


120 


1.6611e-03 


1.9555 


5.7775C-04 


1.9011 


3.6447e-04 


1.6804 




0.2 


240 


4.3643e-04 


1.9283 


1.5215e-04 


1.9250 


1.0389e-04 


1.8107 






480 


1.1223e-04 


1.9593 


3.9170e-05 


1.9577 


2.7629e-05 


1.9109 






60 


7.5867e-03 


- 


2.4101e-03 


- 


1.3364e-03 


- 






120 


2.0069C-03 


1.9185 


6.4998C-04 


1.8906 


3.7650e-04 


1.8277 




1 


240 


5.1832e-04 


1.9531 


1.6801C-04 


1.9519 


1.0062e-04 


1.9037 






480 


1.3136C-04 


1.9803 


4.2497C-05 


1.9831 


2.5599e-05 


1.9748 






60 


1.1959e-02 


- 


3.8026C-03 


- 


1.9938e-03 


- 


Ub 


= a 


120 


3.2940e-03 


1.8602 


9.9527e-04 


1.9338 


5.4231e-04 


1.8783 


T = 


5 


240 


8.7736e-04 


1.9086 


2.5358e-04 


1.9727 


1.3933e-04 


1.9606 






480 


2.3271e-04 


1.9146 


6.4252e-05 


1.9806 


3.4967e-05 


1.9944 






60 


5.7714e-03 




1.9358e-03 




1.0481e-03 




Ub 


= 0.75 


120 


1.5035e-03 


1.9406 


5.1617e-04 


1.9070 


2.8061e-04 


1.9011 


T = 


0.2 


240 


3.9299e-04 


1.9357 


1.3616C-04 


1.9225 


7.9134e-05 


1.8262 






480 


1.0063e-04 


1.9655 


3.5080C-05 


1.9566 


2.1035C-05 


1.9115 






60 


7.1823e-03 




2.2843e-03 




1.2069C-03 




Ub 


= 0.75 


120 


1.8963e-03 


1.9213 


6.1315C-04 


1.8974 


3.4013e-03 


1.8272 


T = 


1 


240 


4.8284e-04 


1.9736 


1.5796C-04 


1.9567 


9.0912e-04 


1.9035 






480 


1.2093e-04 


1.9974 


3.9783e-05 


1.9894 


2.3121e-05 


1.9753 






60 


1.1042e-02 




3.5020e-03 




1.8299e-03 




Ub 


= 0.75 


120 


3.0287e-03 


1.8662 


9.1181e-04 


1.9414 


4.8976e-04 


1.9016 


T = 


5 


240 


8.0111e-04 


1.9186 


2.3118e-04 


1.9797 


1.2593e-04 


1.9595 






480 


2.1076e-04 


1.9264 


5.8358e-05 


1.9860 


3.16276-05 


1.9934 



The accuracy test for the 
1 and M = 



Table 4.1 

MBL equation (|4.1[) with e 



trapezoid scheme for the 
= 2. 



352 
353 
354 
355 
356 
357 



|4.2| shows that the semi-discrete scheme has the order of accuracy greater than 2.5 
for ah the cases, and exceeds 3 for some cases. This confirms the accuracy study 
given in sections 3.1.1 and 3.2 respectively. 



We will now use examples to study the solutions to MBL equation (4.1 1 using 
the numerical schemes proposed in chapter [3j We first notice that if we scale t and 
X as follows 

~ t ^ X 



then MBL (4.1) equation can be written in terms of t and x as follows 
(4.3) Ui + = + Tu^^i- 
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AT 


II 

II 2 


order 

1 


II 

^Ax — M Ax 
II 2 


order 

2 


II 

^Ax — M Ax 
II 2 


order 






120 


2.6992e-03 


- 


1.1300C-03 


- 


7.2363e-04 


- 




= 0.9 


240 


4.0403e-04 


2.7400 


1.7079e-04 


2.7260 


1.1283e-04 


2.6811 




0.2 


480 


5.7504e-05 


2.8127 


2.4624C-05 


2.7941 


1.6242C-05 


2.7963 






960 


8.4934e-06 


2.7592 


3.0892e-06 


2.9948 


1.7607e-06 


3.2055 






120 


4.7731e-03 


- 


2.0192e-03 


- 


1.7267e-03 


- 




= 0.9 


240 


8.7205e-04 


2.4524 


3.6879e-04 


2.4529 


3.0632e-04 


2.4949 




1 


480 


1.2006e-04 


2.8606 


5.0480e-05 


2.8690 


4.1985e-05 


2.8671 






960 


1.5942e-05 


2.9129 


6.6663e-06 


2.9208 


5.1464e-06 


3.0282 






120 


3.7573e-03 


- 


1.2122e-03 


- 


7.9211C-04 


- 




= 0.9 


240 


7.4624e-04 


2.3320 


2.4164e-04 


2.3267 


1.5061C-04 


2.3949 




5 


480 


1.1994e-04 


2.6373 


3.8434e-05 


2.6524 


2.5089e-05 


2.5857 






960 


1.5565e-05 


2.9460 


4.9190e-06 


2.9660 


3.1363e-06 


2.9999 






120 


2.1836e-03 


- 


9.1039e-04 


- 


5.7219e-04 


- 




— fj 


240 


3.2729e-04 


2.7381 


1.3760e-04 


2.7260 


8.9550e-05 


2.6757 




0.2 


480 


4.6856e-05 


2.8043 


1.9909e-05 


2.7890 


1.2935e-05 


2.7914 






960 


6.7382C-06 


2.7978 


2.3182e-06 


3.1023 


1.4109e-06 


3.1965 






120 


3.9014C-03 


- 


1.6388e-03 


- 


1.3873e-03 


- 






240 


7.0517C-04 


2.4680 


2.9669e-04 


2.4656 


2.4272e-04 


2.5149 




1 


480 


9.6528e-05 


2.8690 


4.0354e-05 


2.8781 


3.3125e-05 


2.8733 






960 


1.2890e-05 


2.9047 


5.3648e-06 


2.9111 


4.07546-06 


3.0229 






120 


3.0797e-03 


- 


9.9202e-04 


- 


6.4456e-04 


- 


Ub 


= a 


240 


6.1133C-04 


2.3328 


1.9783C-04 


2.3261 


1.2277e-04 


2.3924 


T = 


5 


480 


9.7351e-05 


2.6507 


3.1222e-05 


2.6637 


2.02636-05 


2.5990 






960 


1.2396e-05 


2.9733 


3.9513e-06 


2.9822 


2.49626-06 


3.0210 






120 


1.8244e-03 




7.5548e-04 




4.66716-04 




Ub 


= 0.75 


240 


2.7262e-04 


2.7425 


1.1419e-04 


2.7260 


7.32996-05 


2.6707 


T = 


0.2 


480 


3.9198e-05 


2.7980 


1.6562e-05 


2.7855 


1.06816-05 


2.7788 






960 


5.4739e-06 


2.8401 


1.9677e-06 


3.0733 


1.3232e-06 


3.0129 






120 


3.2727e-03 




1.3672C-03 




1.1477e-03 




Ub 


= 0.75 


240 


5.8671e-04 


2.4798 


2.4585e-04 


2.4754 


1.98666-04 


2.5304 


T = 


1 


480 


7.9974e-05 


2.8750 


3.3285e-05 


2.8848 


2.70336-05 


2.8775 






960 


1.0724e-05 


2.8987 


4.4466e-06 


2.9041 


3.33416-06 


3.0193 






120 


2.5902e-03 




8.3335e-04 




5.38826-04 




Ub 


= 0.75 


240 


5.1342e-04 


2.3348 


1.6611e-04 


2.3268 


1.02716-04 


2.3913 


T = 


5 


480 


8.1062e-05 


2.6630 


2.6032e-05 


2.6738 


1.68136-05 


2.6109 






960 


1.0173e-05 


2.9944 


3.2662e-06 


2.9946 


2.04736-06 


3.0377 



Table 
scheme 



for the MBL equation ([4J 



for the third order semi-discrete 
11) with e = 1 and M = 2. 



359 
360 
361 



364 
365 



The scaled equation (4.3 1 shows that it is the magnitude of ^ and 

X 



that determine 



the asymptotic behavior, not t, x, neither e alone (|2T]). In addition, ( 4.3 1 also shows 
that the dispersive parameter r denotes the relative importance of the dispersive 
term u^xt- The bigger r is, the more dispersive effect (4.1) equation has. This can 



363 be seen from the computational results to be shown later in this section. 



Duijn et al. [3T] numerically provided a bifurcation diagram (Figure 4.2 1 of MBL 



(4.1) equation as the dispersive parameter r and the post-shock value ub of the 

366 initial condition vary. The solution of (4.1 1 has been proven to display qualitatively 

367 different profiles for parameter values {t,ub) falling in different regimes of the 

368 bifurcation diagram. In particular, for every fixed r value, there are two critical 
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bifurcation diagram 




Figure 4.2. The bifurcation diagram of the MBL equation (1.10 1 
with the bifurcation parameters {t,ub)- 



369 
370 
371 

372 
373 
374 

375 
376 



379 
380 
381 
382 
383 



389 
390 
391 
392 
393 
394 



Ub values, namely, u and u. From the bifurcation diagram (Figure 4.2), it is clear 
that, when t < t^,, u = u ~ a. For a fixed r value, the solution has three different 
profiles. 

(a) If Ub G [u, 1], the solution contains a plateau value ub for < | < ^("s), 
a rarefaction wave connection ub io u for ^{ub) < f < ^(^)^ another 
plateau value u for ^{u) < | < Ar^, and a shock from u down to at 
f — (see Figure 

(b) If ub G (u, u), the 
•^^"j"^'-""'' , a shock from ub up to S at 
value n for < ? < ^ 

X _ f{u) 



3(a)) 



solution contains a plateau value for < f < 

another plateau 



f(u)-f(uB) 



t u—ub 

and a shock from u down to at 



(see Figure 3(b) I. The solution may exhibit a damped oscillation 
near u — ub- 

If ub G (0,u], the soluti on co nsists a single shock connecting ub and 



at 



(see Figure 3(c)). It may exhibit oscillatory behavior near 



u — ub- 

384 Notice that when r > and u < ub < u, the solution profiles ( |3(b)[ ) displays 

385 non-monotonicity, which is consistent with the experimental observations ([7]). 

386 In the numerical computation we show below, we will therefore test the accuracy 

387 and capability of central schemes for different parameter values (r and ub) that fall 

388 into various regimes of the bifurcation diagram, and therefore display qualitatively 
different solution profiles. The numerical experiments were carried out for M = 2, 
e = 0.001 and T = 4000 x e, i.e. T = 4000 to get the asymptotic solution profiles. 



and Ax was chosen to be and A 



^ was chosen to be 0.1. The scheme used 



in the computation is the second order Trapezoid scheme as shown in section [3. 1.1| 
The Midpoint scheme delivers similar computational results, hence is omitted here. 



The solution profiles at j (blue). 



2*T 
4 



(green), 



3*T 



(magenta) and T (black) are 
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(a) 



(b) 



(c) 



Figure 4.3. Given a fixed r, the three quahtatively different so- 
lution profiles due to different values of Ms. In particular, when 
T > r* and u < ub < u, the solution profiles (Figure |3(b)| ) dis- 
plays non-monotonicity, which is consistent with the experimental 
observations ([7]). Figures [3(a)| [3(b)] and [3(c)] are demonstrative 
figures. 



395 chosen to demonstrate the time evolution of the solutions. The red dashed lines are 

396 used to denote the theoretical shock locations and plateau values for comparison 

397 purpose. 

398 We start with t > 0. Based on the bifurcation diagram (Figure 4.2), we choose 

399 three representative ub values, i.e. ub — 0.9 > a, ub = a = \f^^^ = \/| (for 

400 
401 

402 u < Ub = 0.9, and t = 5 with ub = 0.75, a, 0.9 G [u^^s, Ut-^s]. We first use this 9 

403 pairs of {t,ub) values given in Table 4.3 to validate the solution profiles with the 
demonstrative solution profiles given in Figure |4.3[ 



M+l V 3 

M = 2) and ub = 0.75 < a. For each fixed ub, we choose three representative r 
values, i.e. r = 0.2 < « 0.61, r = 1 > r* with ub = 0.75 < Ut^i < ub = a < 



(t, Ub) 


Example 4 


Example 5 


Example 6 


Example 1 


(0.2,0.9) 


(1,0.9) 


(5,0.9) 


Example 2 


(0.2, a) 


(1,«) 


(5, a) 


Example 3 


(0.2,0.75) 


(1,0.75) 


(5,0.75) 



Table 4.3. 9 pairs of {t,ub) values with either fixed t value or 
fixed Ub value used in Examples 1-6. 



407 
408 
409 
410 



Example 1 (t, us) ^ (0.2, 0.9), (r, ub) = (1, 0.9), (r, ub) = (5, 0.9). 
When ub = 0.9 > a is fixed, we increase r from 0.2 to 1 to 5 (Figure [4(a)[ , |4(b)| 



4(c) I, the dispersive effect starts to dominate the solution profile. When r — 0.2 



(Figure 4(a)), the solution profile is similar to the classical BL equation solution 
(see Figure 2(b) ), with a rarefaction wave for j E [f'{u = 0.9), f'{u — a) ^ f'{u = 
Ut=o.2)] and a shock from u = cktoM = Oat| = /'(o^). This corresponds to 



3(a) with ^(Mr=o.2 = a) 



= ^^^;=°f = When T = 1 (Figure |4(b)| , 

S [f'{u = 0.9), /'(u = Ut=i)] and the solution 

^^^^1 and the shock 



411 Figure 

412 the rarefaction wave is between 

413 remains at the plateau value u — Ut=i for | G [/'(« = Ur=i), 
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414 

415 
416 

417 
418 



occurs at 



_ /(«T=1 



This corresponds to Figure 3(a 



with ub — 0.9 > Ur=i 



0.86. When r = 5 (Figure [4(c) ), the solution displays the first shock from u — 0.9 

X _ /(«t = 5)-/(Mb) 

r /(«x=5)-/(«fl) /(°^=5) j ^j^^ ^j^g second shocks occurs at f 



to u = Ut-^s at 
for f e 



and then remains at the plateau value u — Ur=5 



/("t = 5) 



This 



corresponds to Figure [3(b) with Ut=5 ~ 0-68 < ub = 0.9 < Mr=5 ~ 0.98. Notice 
419 that as T increases, the rarefaction region shrinks and the plateau region enlarges. 



(a) {r,«s) = (0.2,0.9) 
1 



=> 0.5 



(d) (t,«s) = (0.2, a) 



=1 0.5 



2 4 

X 

(g) (T,«fl) = (0.2,0.75) 
1 

I I I 



0.5 



(b) (t,«s) = (1,0.9) 



=1 0.5 



2 4 

X 

(e) {t,ub) = {I, a) 



=3 0.5 



2 4 

X 



(h) (t,ub) = (1,0.75) 



=3 0.5 




(c) {t,ub) = (5,0.9) 



=i 0.5 



T~TTT 



(f) (r,«s) = (5, a) 



=1 0.5 



2 4 

X 

(i) {t,ub) = (5,0.75) 
1 



=3 0.5 



Figure 4.4. Numerical solutions to MBL equation with param- 
eter set ting s fall in different regimes of the bifurcation diagram 
(Figure 4.2 1. The color coding is for different time: jT (blue), |r 



(green), |T (magenta) and T (black). The results are discussed in 

examples 1-6. In figures 
M = 2. 



4(d) 



4(f) 



= for 

M+l V 3 



421 Example 2 (r, ub) = (0.2, a), (r, ub) = (1, a), {t,ub) — (5, a). 

422 When Ub — a is fixed, we increase r from 0.2 to 1 to 5 (Figure [4(d)| , |4(e)| , |4(f)[ ), 

423 the dispersive effect starts to dominate the solution profile. When r = 0.2, the 



424 solution displays one single shock at 



For both r = 1 and r = 5, the 
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427 
428 
429 
430 
431 
432 



439 
440 



448 
449 
450 
451 
452 
453 



459 
460 
461 
462 
463 



469 
470 

471 



solution has two shocks, one at | 



/("t=1(t = 5 



oly))-/(") 



426 at 



/(^T = l(-r = 5 respectively)) 

'^7- = l(-r = 5 respectively) 



"t=1(t = 5 respectively) —a 

For both t — 1 and t — 5 (Figures 



and another one 
the 



4(e) 



4(f)) 



solutions correspond to Figure [3(b)| which are consistent to the experimental obser- 
vations. Notice that as r increases from 1 to 5, i.e., the dispersive effect increases, 
the inter-shock interval length increases at every fixed time (compare Figure |4(e) 
with Figure [4(f)] ). In addition, for fix t = 1 (t = 5 respectively), as time progresses, 
the inter-shock interval length increases in the linear fashion (see Figure 4(e) (Fig- 
ure [4(f)] respectively) ). 

Example 3 (r, ub) = (0.2, 0.75), (r, us) = (1, 0.75), (t, ub) = (5, 0.75). 
When Ub = 0.75 <= a is fixed, we increase r from 0.2 to 1 to 5 (Figure 4(g)[ , [4(h)] 



4(i)|, the dispersive effects starts to dominate the solution profile in the similar 
0.9 and ub — a. Notice that when t — 1, since ub 

(Figure 4(h)[) 



fashion as ub 
very close to Ur=i, 



a. Notice that when r 
the solution displays oscillation at 



0.75 is 
If 



we increase r further to r = 5, the dispersive effect is strong enough to create a 
plateau value at u 



0.98 (see Figure 4(i) ) 



Example 4 (t,ub) = (0.2, 0.9), (r, wb) = (0.2, a), (r, ws) = (0.2 0.75). 

Now, we fix T = 0.2, decrease ub from 0.9 to a, to 0.75 (Figure^ [4(a)]|4(d)[[4(g)| ). If 
Ub > ot the solution consists a rarefaction wave connecting ub down to a, then a 
shock from a to 0, otherwise, the solution consists a single shock from ub down to 
0. In all cases, since r = 0.2 < r*, regardless of the ub value, the solution will not 
display non-monotone behavior, due to the lack of dispersive effect. 

Example 5 (r, u^) = (1, 0.9), (r, ub) = (1, a), (t, ub) = (1, 0.75). 
Now, we fix r = 1, decrease ub from 0.9 to a, to 0.75 (Figure^ [4(b)][4(e)][4(h)] ) . If 
Ub — 0.9 > Ut=1i the solution consists a rarefaction wave connecting ub and u, 
and a shock connecting u down to (Figure [4(b)] ). Even if u < Ub < u, because 
T = 1 > , the solution still has a chance to increase to the plateau value u as seen 
But, if Ub is too small, for example, ub 



454 in Figure 4(e) 



0.75 < u, the solution 
does not increase to u any more, instead, it consists a single shock connecting ub 
down to (Figure [4(h)[ ) . 

Example 6 (r, us) = (5, 0.9), (r, ub) = (5, a), (r, ub) = (5, 0.75). 

Now, we fix T = 5, decrease ub from 0.9 to a, to 0.75 (Figures 4(c)|4(f)|4(i)| ). For all 
three ub, they are between Mt=5 ^-nd Ur=5, hence all increase to the plateau value 
Ur=5 sa 0.98 before dropping to 0. Notice that as ub decreases, the inter-shock 
interval length decreases at every fixed time (compare Figures [4(c)[ [4(f)] and 4(i) ). 
This shows that when the dispersive effect is strong (t > r*), the bigger ub is, the 
bigger region the solution stays at the plateau value. 



Example 7 (r, ub) = (0, 0.9), (r, ub) = (0, a), (r, ub 
We now show the solution profiles for the extreme r value, i.e 
0.9), [5(b)] (hb 



5(a) {ub = 0.9), [5(b)] {ub = a) and 5(c) {ub = 
of classical BL equation with small diffusion eu 
and 



0.75). 



(0,0.75). 

r = in Figures 
Notice that these are cases 



2(b) 



5(c) with the solution of the classical BL equation given in Figures|2(a) 



5(a)| 


5(b) 


s|2(a) 


and 



it is clear that they show qualitatively same solution profiles. The difference 
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472 

473 



is that due to the diffusion term in the MBL equation, as shown in Figure [475] the 
sohitions do not have sharp edges right at the shock, instead, the solutions smear 
out a httle. If we compare Figures [5(a)| [5(b)] and [5(c)] with Figures [4(a)] [4(d)| and 



4(g) there is no visible difference. This shows that once r < t«, solution profile 



will stay the same for a fixed ub value. 



(a) {t,ub) = (0,0.9) 




(b) (r,ns) = (0,a) 



=> 0.5 




(c) {t,ub) = (0,0.75) 



=> 0.5 



0.5 



Figure 4.5. The numerical solutions of the MBL equation at T 
= 1 with T — and different ub values. The results are discussed 
in example 7. 
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Example 8 (t,ub) = (5, 0.99), (r, Uij) = (5, 0.98), (t, us) = (5,0.97). 
We also study the solution profiles for ub close to u. For example, when r = 5, 
u « 0.98, we hence choose ub = 0.99, ub = 0.98, ub = 0.97 and solutions are 
shown in Figure 6(a)[ |6(b)[ [6(c) If ub = 0.99 > Ur=5 ~ 0.98, the solution drops 
to the plateau value u, then drops to (see Figure 6(a)[ ). If ub — 0.98 w Mr=5, 
the solution remains at plateau value Ur=5 and then drop to (see Figure [6(b)| . 
If Ub = 0.97 < Ut-=5, the solution increases to the plateau value Ur=5 ~ 0.98, then 
drops to 0. In all cases, the transition from ub to Ut=5 ~ 0.98 takes very small 
space. In the majority space, the solution keeps to be the plateau value Ur=5 ~ 0.98. 



(a) {t,ub) = (5,0.99) 
1 r ■ ■ — 



0.99 
0.98 
0.97 
0.96 



0.2 0.4 O.f 

X 



(b) (t,ub) = (5,0.98) 

1 r ■ ■ — 

0.99 I 
=> 0.98 1 



0.97 
0.96 



0.2 0.4 0.6 

X 



(c) {t,ub) = (5,0.97) 
1 r ■ ■ — 



0.99 
=3 0.98 
0.97 
0.96 



0.2 0.4 0.6 

X 



Figure 4.6. Numerical solutions to MBL equation with ub close 
to Ur=5 « 0.98. The color coding is for different time: |T (blue), 
|r (green), |T (magenta) and T (black). The results are discussed 
in example 8. 



488 Example 9 {t,ub) ~ (5,0.7), {t,ub) — (5,0.69), {t,ub) = (5,0.68), {t,ub) 

489 (5,0.67), (r,UB) = (5,0.66). 
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491 



496 
497 
498 



In addition, we study the solution profiles for ub close to u. For example, when 
T = 5, u sa 0.68, we hence choose ub — 0.7, ub — 0.69, ub = 0.68, ub = 0.67, 
Ub = 0.66 and solutions are shown in Figures 7(a)[ |7(b)[ |7(c)[ |7(d)[ |7(e) As 
Ub decreases crossing u^^c, ~ 0.68, the solution gradually stops increasing to the 
plateau value Ut=5 , and the inter-shock interval length decreases (compare Figures 



7(a)[ |7(b)| and |7(c)[ ). The oscillation in Figures |7(d)| and [7(e) are due to the fact 
that Ub values are too close to u^.^^- This confirms that even with big dispersive 
effect (say r = 5), if Ms is too small (e.g. ub < u), the solution will not exhibit 
non-monotone behavior. 



(a) (t,«s) = (5,0.7) 



0.5 



(b) (r,ns) = (5,0.69) 



=1 0.5 



(c) {t,ub) = (5,0.68) 



=3 0.5 



(d) (t,«s) = (5,0.67) 



(c) (r,«fl) = (5,0.66) 



=1 0.5 



=! 0.5 



Figure 4.7. Numerical solutions to MBL equation with ub close 
to Ut=5 ~ 0.68. The color coding is for different time: |T (blue), 
|T (green), |T (magenta) and T (black). The results are discussed 
in example 9. 
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Example 10 {t,ub) = (0.2,0.6), {t,ub) = (1,0.6), {t,ub) = (5,0.6). 
We fix Ms to be small, and in this example, we take it to be ub = 0.6. We vary the r 
value, from t = 0.2 < t^, to t — 1 barely larger than to r = 5 > . The numerical 
solutions are given in Figure |8(a)[ |8(b)[ |8(c) As t increases, the post-shock value 
remains the same, but there will be oscillation generated as r becomes larger than 
. Figures [8(d)| [8(e)] and [8(f)] show that as r increases, the oscillation amplitude 
increases and oscillates more rounds. Notice that r is the dispersive parameter, and 
this means that even for small ub value, different dispersive parameter values still 
give different dispersive effects, although none can bring the solution to the plateau 
value u. Comparing Figures [8(d)] [8(e)] and [8(f)] with Figures [8(g)] [8(h)] and [8^ 
is clear that the oscillation amplitude remains steady with respect to time. 
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513 



Example 11 e = 0.001, e = 0.002, e = 0.003, e = 0.004, e = 0.005. 

In this example, we will compare the solution profiles for different e values. Fixing 
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(a) (r,«s) = (0.2,0.6) 
1 



0.5 



I I 
I I 



I 



(b) {t,ub) = (1,0.6) 



=1 0.5 



(c) (t,«b) = (5,0.( 



=> 0.5 



(d) Fig 8(a) zoomed in at (e) Fig 8(b) zoomed in at |T (f) Fig 8(c) zoomed in at |T 



=1 0.5 




=1 0.5 



0.84 0.86 0.88 0.9 0.92 

X 




=1 0.5 




0.84 0.86 0.88 0.9 0.92 

X 



0.84 0.86 0.88 0.9 0.92 

X 



(g) Fig |8 (a) I zoomed in at T (h) Fig|8(b)|zoomed in at T (i) Fig|8(c)|zoomed in at T 



=1 0.5 




=1 0.5 



3.48 3.5 3.52 3.54 3.56 




=1 0.5 



3.48 3.5 3.52 3.54 3.56 




3.48 3.5 3.52 3.54 3.56 



Figure 4.8. Numerical solutions to MBL equation with small con- 
stant ub — 0.6 and different t values. The figures on the second 
and third rows are the magnified versions of the first row at t = 
and t = T respectively. The color coding is for different time: \T 
(blue), |T (green), |T (magenta) and T (black). The results are 
discussed in examples 10. 



T = 0.5, Aa; ~ 0.0001, A = ^ = 0.1, we show the numerical results in Figure 



4.9 



615 for e = 0.001 (blue), e = 0.002 (yellow), e = 0.003 (magenta), e = 0.004 (greell); 

516 and e = 0.005 (black). For the purpose of cross reference, we choose the same 

517 nine sets of parameter settings as in examples 1- 6. To assist the observation, the 

518 figures in Figure |4.9| are zoomed into the regions where different e values introduce 

519 different solution profiles. The numerical solutions clearly show that as e increases, 

520 the numerical solution is smeared out, and the jump location becomes less accurate. 

521 Notice that r is responsible for the competition between the diffusion and disper- 

522 sion, which in turn determines the plateau values. Hence varying e value doesn't 

523 affect the plateau location. 
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(a) (r,«s) = (0.2,0.9) 
1 I 




(d) (t,«s) = (0.2, a) 




(g) (r,UB) = (0.2,0.75) 
1 



=> 0.5 




(b) (r,«s) = (1,0.9) 



=1 0.5 




0.5 



0.55 

X 



0.6 



(c) {t,ub) = {I, a) 



=1 0.5 



0.4 0.5 0.6 

X 

(h) {t,ub) = (1,0.75) 



=1 0.5 




(c) (t,«b) = (5,0.9) 



=> 0.5 



0.2 0.4 0.6 

X 

(f) (r,«B) = (5,a) 



=> 0.5 



0.2 0.4 0.6 

X 

(i) {t,ub) = (5,0.75) 



=> 0.5 




Figure 4.9. The numerical solutions of MBL equation at T = 0.5 
with e = 0.001 (blue), e = 0.002 (yellow), e = 0.003 (magenta), 
e = 0.004 (green), and e — 0.005 (black). The view windows are 
zoomed into the regions where different e values impose different 
solution profiles. The results are discussed in example If. 
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5. Conclusion 

We proved that the solution to the infinite domain problem can be approximated 
by that of the bounded domain problem. This provides a theoretical justification 
for using finite domain to calculation the numerical solution of the MBL equation 
(f.lO). We also extended the classical central scheme originally designed for the 



hyperbolic systems to solve the MBL equation, which is of pseudo-parabolic type. 
The numerical solutions for qualitatively different parameter values r and initial 
conditions ub show that the jump locations are consistent with the theoretical 
calculation and the plateau heights are consistent with the numerically obtained 
values given in |2f j . In particular, when t > r*, for ub G (m, u), the numerical 
solutions give non-monotone water saturation profiles, which is consistent with the 
experimental observations. In addition, the order tests show that the proposed 
second and third order central schemes achieved the desired accuracies. 
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538 In [221 EQ], the two-dimensional space extension of the modified Buckley-Leverett 

539 equation has been derived. One of the future directions is to develop high order 

540 numerical schemes to solve the two-dimensional MBL equation. Central schemes 

541 have been used to solve high dimensional hyperbolic problem and dispersive prob- 

542 lem ([Ullin]), which makes it a good candidate for such a task. 
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544 



Appendix A. Proof of the lemmas 



Proof to lemma\2^ Let g{u) = ^ = ^2j^m"{i-uY 



then 



, _ M -{1 + M)u^ 



>0 if Q<u<J^^ 



= if u = 



M 
M+l 



<0 if u>J^^ 



545 and hence g(u) achieves its maximum at u = y mXi • Therefore, = 5; D, 

546 where D — ^-^^ and a = y j^^, and in turn, we have that f{u) < Du for all 



547 < It < 1. 

Proof to lemma Jlp 



□ 



+00 



e — e 



e = e^/T- 



-2 + 2e ^ 2e\/r 



A2 - 1 - 1 - A2 



if A e (0,1) 



□ 



Proof to lemma 2.3 ^ 







e — e 



e «^ d(, — xe < 



e(l - A) 



if A e (0,1). 



□ 



550 Proof to lemma 2.3 {Hi). Based on the assumption on uq in (2.27) 

f -t-oo 



(A.l) 



+ 00 



e '^^/^'e«'^^|^io(e)|dC 



551 Calculating yi{x) with the assumption that A G (0, 1), we get 



^ Jo e^v/^ < ey/re < e^ 



for X £ [0, Lq] 
for X G [Lq, +oo) 



552 Therefore, we get the desired inequality 

f +00 







e — e 



lMo(?)Me<2af 



□ 



BOUNDED DOMAIN PROBLEM FOR THE MODIFIED BUCKLEY-LEVERETT EQUATION33 



Proof to lemma\KM fl). 



+00 





A2 - 1 



e + sgn(a; — ^)e 



e 



(A-1)^ 

-2 + 2Ae '=v^ 



2(A- l)e" 



< 



1-A2 



if A e (0,1) 



557 Proof to lemma 2.4 

558 

+ C30 



+ sgn(a; — ^)e 



(A-3)x {X-l)x 

2e - 2e 



^ (A-l)x 



-2 



xe < evr - 



e(l-A) 



if A e (0,1). 



560 Proof to lemma 2.4 {in). Based on the assumption on uq in (2.27) 



+00 



e + sgn(x — ^)e e=^|Mo(C)|c?^ 



(A.2) 



e + sgn(a; — ^)e 



< Cue^ 

561 Calculating 2/3(0;) with the assumption that A e (0, 1), we get for x € [0, Lq] 

(A-l)x J_ ; (A + l)x /""U {_ ALn 



2/3(3^) < 
562 and 



' L (A + l)x 

(e ■'v^ + e-'v^) + e -v/^ 



(a-i)x+j:o 



2/3(2;) < e / (e + e^^) d£, < cy/re < e-\/re=^ 



563 for X e [io, +00). 

564 Therefore, we get the desired inequality 

'■ + °° _x + 5 x-;| Ax , , ^, ^ Aijj 

e + sgn(a; - C)e e--^ \uo{^)\d^ < 2Cue^/Te'-^ 





Proof to lemma \275\ 

4>i{x) — e 



e — e ' 

L 

f ^ - — e ' 



567 Proof to lemma 2.5 Q). Since 02(x) 

568 and hence 02 (a;) < <j>2{L) = 1 for a; G [0,L]. 



, we see that 02 (x) 



□ 



□ 



□ 



□ 

1 e +e 



□ 
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Proof to lemma 2.5 (Hi). (j)2{x) 



■ +e 



gives that (1)2 (x) = -^cj)2{x) > 



570 0, and hence 02 (a;) < 02 (i) = 

671 X e [0, L]. 



• +e 



1 e< 



< if e < 1 for 



□ 
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